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The characteristic formalism in numerical relativity, which has been developed to study grav- 
itational waves, and the observer metric approach in observational cosmology both make use of 
coordinate systems based on null cones. In this paper, these coordinate systems are compared and 
it is then demonstrated how characteristic numerical relativity can be used to investigate problems in 
observational cosmology. In a numerical experiment using the characteristic formalism, it is shown 
how the historical evolution of a LTB universe compares to that of the ACDM model given identical 
observational data on a local observer's past null cone. It is demonstrated that, at an earlier epoch 
of the LTB model, the observational data would not be consistent with that of the ACDM model. 
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I. INTRODUCTION 



Provided that General Relativity is valid up to the largest scales, under ideal circumstances, 
the space-time structure of the Universe can be determined by solving the full set of Einstein field 

q-i equations (EFEs) using observational data as the boundary conditions. Such an observational ap- 

proach makes use of minimal assumption and observations as the boundary conditions dictate the 

(2J[)| outcome of the solution. Because of difficulties in obtaining and interpreting observational data, 

the conventional approach to cosmology is rather based on parameterized models where a model 
is assumed and then validated against different sets of observations and constraints. Assuming 
a homogeneous universe with a non-zero cosmological constant and with suitable adjustments of 
parameters, predictions using the ACDM model show remarkable correlation with observed cos- 
mological parameters. The assumption of homogeneity, which is an integral part of this approach, 
has the implication that local results from simple models can easily be extrapolated to the whole 
CO J Universe. 

However, doing investigations such as quantifying homogeneity on different scales, testing the 
verifiability of cosmology [1], validating the Copernican principle [2] and determining the metric 
of the Universe Q , a more suitable methodology is one where the EFEs are solved using a general 
metric and boundary conditions derived from observations. Since homogeneity is not assumed, 
conclusions following from the observational approach are necessarily limited to the causally con- 
nected region in the interior of our current past null cone on which observations are located. Since 
the EFEs are solved in a general form, or at least more general than the ACDM model, there are 
more degrees of freedom and also more required boundary conditions. Sufficiently accurate and 
complete observations do not currently exist, however, it is anticipated that the next generation 
of astronomical surveys will begin to provide enough detail to perform observational solutions in 
spherical symmetry. 

Much of the development on observational cosmology was inspired by the seminal paper of 
Kristian & Sachs 4] published in 1966. In their work, observational data was used to determine a 
general metric but since they made use of series expansions, it was necessarily restricted to a region 
close to the observer. The ideas introduced by Kristian & Sachs were the starting point of further 
developments by Ellis, Stoeger, Maartens and others in their Observational Cosmology programme 
with its initial publication in 1985 [l|. In these developments, observational coordinates, based on 
the concept introduced by Temple in 1938 [5j, were implemented to extend the region investigated 
by Kristian & Sachs to higher redshifts. Observational coordinates are based on null geodesies 
on the past null cone (PNC) and are therefore the natural framework on which electromagnetic 
radiation reaches an observer. Solving the EFEs for cosmology then consists of two problems: 
firstly, astronomical observations are used to determine the metric on the local PNC and secondly 



these form the final values of a characteristic final value (CFVP) problem which determines the 
historical evolution of the region causally connected to the PNC (i.e. the interior of the PNC). 
The causally connected region is of fundamental importance since it defines the limits on which 
cosmological models can be validated from direct observations. 

In developments towards exact solutions in observational cosmology, spherical symmetrical so- 
lutions are taken as a first step to refine the methods (sec for instance _Q-8j). Spherical symmetry 
by itself is not necessarily unrealistic since the Universe does appear to be highly isotropic. The 
spherical symmetrical inhomogeneous Lemaitre-Tolman-Bondi (LTB) model in observational co- 
ordinates is therefore an important tool for verifying these developments. The LTB model in its 
standard form also provides a useful framework to investigate solutions from direct observations. 
In this approach, observations on the PNC are transformed to cosmological coordinates and then 
related to the coefficients of the LTB metric (e.g. see [9J]). These transformations require numeri- 
cal solutions to handle observational data and recent work by Lu & Hellaby [1(| and McClure & 
Hellaby [ll| developed and refined methods for setting up the local PNC from realistic observations 
with the intention to be implemented on data obtainable in the next generation of astronomical 
surveys. It should be noted that the emphasis in observational approaches lies in the fact that the 
initial (final) conditions are provided by direct observations on the PNC and that the EFEs are 
implemented in a general form with minimal assumptions. 

In numerical relativity (NR) methods to implement general solutions of the EFEs have been well 
established for strong gravitational scenarios where gravitational waves are expected to be formed. 
Similar to electromagnetic radiation, gravitational radiation also propagates along null geodesies 
and null cones also provide a natural frame of reference. In this case, however, it is the future null 
cone that is of interest in the form of a characteristic initial value problem (CIVP). The character- 
istic formalism in NR is based on the theoretical developments of Bondi, van der Burg & Metzer 
[12| and Sachs [131 ] which were part of the Gravitational Waves in General Relativity programme 
initiated by Bondi in the late 1950s. These developments were of fundamental importance to the 
understanding of gravitational waves and with the advancement of computational technology, it 
was recognised that the characteristic formalism holds several computational benefits. Among 
these are: the fact that the EFEs simplify to ordinary differential equations along characteristics, 
which are less expensive to compute and the conformal method, developed by Penrose |14| . can 
be used to represent infinity on a finite grid, making modelling asymptotical behaviour possible in 
full non-linearity. One drawback of a null cone coordinate system is its behaviour around caustics 
where it can become multi- valued and singular. In astrophysical problems, this problem is usually 
avoided by combining the characteristic formalism with the Cauchy formalism where the latter is 
used to solve regions where caustics are expected while the former is used to extract the solution in 
the far field [15j, [l6| . Treating caustics directly on the null cone has also been investigated [17[ but 
this approach has not been implemented numerically. A comprehensive overview of characteristic 
numerical relativity can be found in [18j . 

In the context of cosmology, the CIVP in NR poses the advantages that cosmological develop- 
ments can leverage on previous developments in astrophysical problems, the causally connected 
region of the PNC can be computed using realistic boundary conditions free of constraints and 
the coordinates are simpler than observational coordinates. The simplicity of the coordinates does 
however introduce restrictions in that the location where the PNC refocusses, the apparent horizon 
(AH), becomes multi- valued. The importance of the AH has been emphasised in |19| and Q and 
methods to cross this region have been presented in [H| , [Tljj and [20] and similar approaches will 
probably work for the characteristic formalism as well. Although recognising the importance of the 
AH, the initial steps taken in the current work will be limited to the region where the character- 
istic formalism is well behaved while the region where caustics occur will be postponed for future 
research. Applying the CIVP in NR to observational cosmology has previously been investigated 
by Bishop & Haines in [2lj and the work presented in this paper is a continuation of their work. 
Recent work by Hellaby & Alfedeel |3| defined an algorithm to implement an approach based on 
observational coordinates with the intention of a numerical implementation. In their work, consid- 
erations were taken into account for passing the AH and some other details not taken into account 



in this paper but a numerical implementation of their work has not been presented yet. 

The work presented in this paper is about evolution off the null cone i.e. the second part of the 
observational cosmology problem. We develop a computer code, the input to which is data on the 
PNC obtainable, in principle, directly from observations. The code uses the EFEs to evolve the 
model into the past. The code reported here is a first step towards the construction of a general 
code for observational cosmology, and is limited to the case of spherical symmetry 

The code is used to investigate an issue in cosmology. The standard model of the universe is 
ACDM, but it is well known that the observational data also fits a LTB model. The code is used 
with data on the PNC generated by the ACDM model, and we evolve into the past using the A = 
EFEs. The result is that we compute the data that would be observed on the PNC at an earlier 
epoch within the context of the LTB model. We then ask the question: Could this data also be 
interpreted as being that of some ACDM model? The answer is "No", which indicates that, if the 
universe is LTB, then not only are we at a special position in space, but also that the the present 
time is special in the history of the universe. 

In section |TT1 the CIVP is formulated and related to cosmology. Section HTT1 presents the imple- 
mentation of the code, and section IIVI describes its verification. The results of running the code 
with ACDM model data are given in section |Vj The paper ends with a Conclusion, section I VII 



II. CHARACTERISTIC FORMALISM 

A. Metric and coordinates 

The essence of the characteristic formalism is a frame of reference based on outgoing null cones 
that evolve from values on an initial null cone. The idea is conceptualised in figure [TJ G is a 
timelike geodesic, and u is the proper time on G. Null geodesies emanating from G have constant 
(u, 9, <f), and near G the angular coordinates 9 and <p have the same meaning as in spherical polar 
coordinates. The coordinate r is the diameter, or area, distance from the cone vertex, which means 
that the surface area of a shell of constant r is <lirr 2 . 

A spherically symmetric null cone metric based on the Bondi-Sachs metric is: : 

ds 2 = -e 2(i ( 1 + — J du 2 - 2e 2fi dudr + r 2 {d9 2 + sin 2 9dip 2 }. (1) 

The hypersurface variables, j3 and W , represent the deviation from a Minkowskian null cone and 
the coordinate system is defined such that j3 and W vanish at the vertex of each null cone i.e.: at 
r = 0, equation (fT]) reduces to a Minkowskian metric. 

Substituting ((T|) into the EFEs, using the form R a ^ = 87r(T a (, — \Tg a b), with the stress-tensor 
for a dust-like fluid (T a b — pv a Vb and T = p), leads to expressions for j3 and W: 

/3, r = 2irrp(v 1 ) 2 (2) 

W, r = e 2fi - 1 - 47re 2 V 2 (3) 

with the initial conditions inherent to the coordinate definition as /3(0) = W(0) — 0. Further, 



1 The notation used here is based on that of \'>'A l and substituting W = V — r will give the original notation of 
Bondi in [2j|. 



substituting the dust stress tensor and ([T]) into the conservation equation, T ab b — 0, it follows that: 



«i,u = — \ (2«iWo - V w (vl) 2 )B r + (V W V! - Uo)vi >r + -(^i) 2 K,,, 
v\ I 2 



(4) 
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An equation for Ufj )U is also obtained from this but making use of the normalisation condition, 
g v a Vb = — 1, a direct expression for vo in terms of Ui, /3 and W can be written as: 



vo = -viV w 



2 1 



(6) 



Without any cosmological considerations, having the values on the initial null cone for p and 
Vi, equations © to (JU) forms a hierarchical system that can be solved in the order @, Q and 
©, then solving equations Q and $5§ evolves the system to the next null cone where the process 
can be repeated until the domain of calculation has been covered. Since these equations are all 
dependent on each other an iterative scheme is required for a numerical solution. 



B. Cosmological coordinates 

In observational cosmology, the observer metric is similar to ([l} but allows for a more general 
treatment of the radial coordinate. Using the notation of [6|, the metric is defined as: 



ds 2 = -A 2 (w, y)dw 2 + 2A(w, y)B(w, y)dwdy + C 2 (w, y){d6 2 + sin 2 0dip 2 }. 



(7) 



With reference to the right hand part of figure [TJ the coordinates x % — {w, y, 9, <p} are defined 
as (see [l| and Q): w, the proper time coordinate on an observer's world line (C), such that 
w = constant hypersurfaces are past light cones. The instance wo is the current time of a local 
observer, defining our PNC. y, is a radial coordinate defined as a null geodesic on the past light 
cone comoving with w. There is some freedom in the in the interpretation of y which includes: 
an affine parameter, the redshift z and the diameter distance C(w,y). As with the characteristic 
formalism, 9 and tp are the spherical inclination and azimuth angles respectively 
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FIG. 1: On the left hand side, the null coordinates of the characteristic formalism and on the right, the 
null cone coordinates of observational coordinates. 



The objective is then to determine A(wo,y), B(wo,y) and C{wo,y) from directly observable 
quantities and then to solve a CFVP to determine the interior of the null cone. For the spherically 



symmetrical dust case observable quantities can be limited to: the redshift z, the apparent lumi- 
nosity I and number counts n. The redshift (z) is measured in terms of distance related to I and 
n is obtained in terms of z. These quantities are then related through source evolution functions 
to the absolute luminosity L and mass per source pt. An important aspect here is that source 
evolution functions should be model independent (9J. 

There are three important differences between (jT|) and ([7]). Firstly, the characteristic formalism 
is more restricted in the radial coordinate r which is chosen to be the diameter distance (r = y = 
C(w,y)). This is not a popular choice in cosmology since r will become multi-valued if the null 
cone starts to refocus at the AH, a situation that can be avoided by making use of y as a distance 
on the null cone in terms of an affine parameter v. However, it has been argued in [19| that the 
position of the AH is an observable property by itself. Knowing this, the model can therefore be 
adapted to compensate for refocussing possibly in a similar approach to that used in |l0| . The 
details of such a development for the null cone will however not be considered in the current paper 
and the calculations will be limited to that region of the PNC interior to the refocussing horizon. A 
second aspect of r is that it is not comoving viz. v 1 ^ 0. The third difference is that the direction 
of the null cone is to the future. Changing the direction of the null cone is straightforward by 
recognising that on a radial null geodesic (ds = d6 = dip = and du ^ 0) from the metric {1}: 

Therefore changing the direction of numerical integration can be used to determine the PNC 
instead of a future null cone. 



C. Cosmological observations 



The observational properties on the initial PNC as derived by Bishop and Haines [2l| will be 
presented in a modified form to separate the redshift z and the diameter distance r$. The redshift 
in terms of the luminosity distance (c^) is taken to be directly observable to a sufficiently accurate 
order in the region where the null cone CIVP is well behaved. The reciprocity theorem is used to 
find a relation between z and r (see |24| ) 

z(di) and di, = (1 + z) 2 tq =>■ z = z(ro). (9) 

The redshift is directly related to the time component of the contravariant velocity: 

l + z = v° (10) 

which can be used to determine the covariant velocity v\ : 

v 1 (u ,r) = -e 2 ^v° = -e 2 ^(l + z). (ff) 

The observed number count (n) is a directly measurable function of z and to be used as input 
to the CIVP once it has been related to the proper density p. Since z(ro) is known, n(z) can be 
rewritten as n(ro) — n(z(ro)). It was shown in }21| that the proper number count can be written 
as: 

N - (T^f -" < 12 > 

In terms of the diameter distance, the proper number count can be written as: 



The proper density is then related to the proper number count: 

p(u ,r ) = f(N). (14) 

The details of this relation will not be considered at this stage but in principle it must take into 
account aspects such as dark matter and source evolution, preferably with factors independent of 
an already assumed cosmological model. 

Solving the CIVP from observational quantities, as with the observer coordinate method, requires 
that the values on the initial null cone be determined from observational quantities. These values 
will then be used as the initial values for the evolution into the local past null cone. The solution 
on the initial null cone follows from equations (fTTj) and (fl~4| but they are dependent on (3 which 
is still unknown. Using equations ©, (jlip and (fl"4"]). a differential equation for j3 which is only 
dependent on observational quantities for the initial null cone can be set up as: 

Ar = 2W(iV)(- e ^(l + z)) 2 . (15) 

This equation can be solved numerically with standard ODE techniques provided that f(N) is well 
behaved. 



III. NUMERICAL IMPLEMENTATION 

A. Code outline 

The code described in this section is based on the general 3D code developed in [22| and [25| but 
significantly simplified for spherical symmetry. The numerical scheme applied is a combination of 
second order finite difference methods based on steps half way between the r and u grid points on 
a regular grid. The objective has been to obtain overall second order convergence in both space 
and time up to a distance reasonably close to the location where refocussing will occur. 

The procedure described here, is to solve equations (J2M]) on a rectangular grid, which is described 
in more detail in section fill Bl Solving the hypersurface equations is done with a central difference 
method on half steps between the r-grid points, using: 

ff3=0J-l + -^Cff?rJ+fl*ri-l) ( 16 ) 

with i being the time step and j the radial step. Here g l r is calculated by substituting known 
values into equations @ and ([3]). In order to solve the evolution equations, (0]) and (JSJ) , their 
general form is notated as: 

ui,-u = F v i and p yU = F p (17) 

and as explicit finite differences on a time half step they are written as: 

fy 1 " 1 = v% + A U F;+ 1/2 and p] +1 = p] + AuF;+ 1/2 . (18) 

Here, n is a time iterator that will approach i. In these equations, the numerical values at the 
point (i,j) are used to evaluate the matter terms while the hypersurface derivatives are obtained 
by substituting the point values directly into the numerical forms of equations ||3J) and (J3J). Ra- 
dial matter derivatives are calculated making use of standard central difference formulae (see for 
instance J2^ p.160-161). 

After setting up a suitable grid (as described in section UlIB I) , the numerical algorithm can be 
summarised in the following steps: 



i. Set the p and v\ initial values on to the initial grid points. These values will, in principle, be 
obtained from observations. 

ii. Calculate /3, W and vq from p and v\ on the initial null cone. 

iii. Calculate Fj the from the values of Vi, p , /3, W and vq. 



iv. Set F. 



n+l/2 



Fj and calculate v i and p as an initial approximation that will approach the 
actual values with subsequent iterations. 

v. Use the new values of v\ and p to calculate /3, W and vq and their radial derivatives, 
vi. Calculate F™ from values in (v) and again v\ and p. 

vii. Test the calculations in vi for accuracy and convergence. If they are sufficiently accurate, 
move to the next time step, otherwise repeat steps v. and vi. with the new values of v\ and 
P- 



B. Grid and domain of calculation 



The domain of calculation is the interior region of the past null cone starting from the present 
space-time location up to a region approaching the point where the null cone starts to refocus. A 
rectangular grid is used to represent the null cone as illustrated in figure [U An important aspect 
of the code is that, besides the coordinate conditions f3(u, 0) = W(u, 0) = 0, the evolution scheme 
is unconstrained on the inner and outer radial boundary i.e. no central or outer conditions will be 
required on the r coordinates. These boundaries do, however, still require to be treated differently 
from the interior region. 
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FIG. 2: The past null cone calculated on a rectangular grid 

On the central world line (r = 0) in figure [5] and its close vicinity, the CIVP equations are 
not suitable in general form (equations (|4|) and (0), which is evident from the occurrence of r 
denominators. Consequently, the r ~ region is calculated by making use of series expansions. In 
order to do this, all the variables in the CIVP equations are replaced by power series expressions 
from which the terms can then be calculated from the resulting coefficients. In the expansions, 
only the terms required to produce second order accuracy will be used, for instance for the density 
evolution: 



P,u = puo + Puir + p U 2r 2 + 0(r 3 ) = F p 



(19) 



The region where the series solution meets with the CIVP solution, also requires special treatment 
to avoid artificial instabilities which has been done by smoothing out the merger region with a 
weighted average between the two solutions. 

In order to evolve the CIVP equations without any boundary condition on the outer radial 
boundary, the outer boundary must be an incoming null hypersurface. This is needed for well- 
posedness, and if violated leads to an numerical instability. An incoming null ray (i.e. on a radial 
null geodesic ds = dO = dip = and du ^ 0) is selected for the outer radial limit and is defined 
from the metric in the following way: 

1 

~ ~2 



?*oute 



'oute 



l 



— ) du 
r 



(20) 



Numerically this can be solved by using an iterative scheme which starts with estimating W = 
and is then repeated until the values of W and r converge. Since the value of r does not necessarily 
fall directly on grid points, W is determined by linearly interpolating between the closest grid points 
on a specific Ui grid line. The behaviour of the density evolution is shown in figure EH where the 
computation extends beyond the null ray cut-off region in order to illustrate the stability issue. 
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FIG. 3: Unstable regions being excluded using a null geodesic on the outer boundary. The figure is based 
on a portion of a normalised Einstein-de Sitter calculation. 



IV. CODE VERIFICATION 



LTB models 



The spherically symmetric CIVP for cosmology is effectively the Lcmaitre-Tolman-Bondi (LTB) 
model in null coordinates. Originating and developed by Lemaitre |27|, Tolman [28| and Bondi [29J . 
the LTB model is a generalisation of the Friedmann-Lemaitre-Robertson- Walker (FLRW) model 
with inhomogeneities allowed in the radial coordinate. The parabolic case of the LTB model is of 
interest for the investigations in this paper and will be used to verify the CIVP code. The metric 
in comoving synchronous coordinates is: 



ds 2 



-dt 2 



[R r (t,r)] 2 dr 2 + [R(t,r)} 2 {de 2 



'} 



(21) 



Here, t is the cosmic (proper) time, r is a comoving radial coordinate with and tp the inclination 
and azimuth angles. R(t, r) is the areal radius and 47ri? 2 defines the proper surface area of a sphere 
with coordinate radius rata constant time slice [3(| ■ Substituting into the Einstein field equations 
and solving gives (see [31|) : 

-l 1/3 



R(t,r) 



-M(r)(t-t B (r)Y 



and p(t, r) 



4tt R 2 R ,, 



(22) 



M(r) is the active gravitational mass which is the mass contributing to the gravitational field and 
is(r) is defined as the bang time function, which is a surface defined by the local time at which 
R = 0. 

In order to provide data to verify the CIVP code against known solutions, LTB models need 
to be transformed to null coordinates which will then be used as input data on the initial null 
cone as well as comparative data inside the null cone. Numerous LTB models exist which include 
the FLRW models, other physical realistic models and also many models with uncertain physical 
relevance. Many physical LTB scenarios have been investigated that are possible candidates of the 
actual universe (see for instance [32|]) which will all be interesting scenarios to run on a CIVP code. 
However, for the purpose of code verification, it will be conducive to work with models that are 
mathematically simple, well documented and physically relevant. The class of LTB models that 
will be investigated is referred to as bang time models and makes use of the bang time function, 
tB, as a mechanism to induce inhomogeneities into a universe. By selecting simple functions for 
t B , it has been demonstrated that inhomogeneities can mimic the effect of inflation (solving the 
horizon problem) [33, |34|i describe the cosmic microwave background (CMB) dipole anisotropy 
|35| . and how inhomogeneities can reproduce the effect of supernovae-redshift dimming without 
the requirement of dark energy [36] . 

For these models, M(r) = Mq = r 3 is chosen as a coordinate condition where Mo is a constant 
which for illustrative purposes is set to 2/9. Equations (|2"2l then reduce to: 



with: 



and p in equation (|22|) to 



R(t,r)=r(t-t B (r)) 2 / ;i (23) 



t-t B (r)-2/3rt B , r {r) , * 

H ^ r >- (t-t B (r))V3 {M) 



P{t ' r) 2n(t - t B (r))(3t - 3t B (r) - 2rt B , r (r)) ' (25) 

By selecting t B (r) = 0, it can be seen that, equation (1251) becomes the Einstein-de Sitter (EdS) 
model, which has been used as a test compared in previous work [2l|. If t B (r) ^ and t B ^ r {r) = 0, 
the time of the initial singularity is adjusted and the age of the Universe changes. For non-constant 
functions, the initial singularity becomes a singular surface and the age of the Universe becomes 
subject to the position of an observer (i.e. the age of the Universe depends on r). Thus, a variety 
of models can be generated for testing the code on parabolic spatial sections. 



B. Coordinate transformations 

In order to generate results from LTB models on the null cone, a coordinate transformation must 
be made from equation ([TJ to equation (|21[) . For simplicity, the null cone metric will be written in 
the form: 

ds 2 = -h u du 2 - 2h r dudr + r 2 {d6 2 + sin 2 6dp 2 }. (26) 

and the coefficients will be rewritten in terms of W and (3 afterwards. Also, LTB variables that 
have the same symbols as the null cone variables will be notated using a tilde e.g. tltb = f, 
vi ltb = vi, g L TB = g etc. 

The essence of the transformations is to relate the (u, r) coordinates to corresponding (t, f) 
values. An expression for r in terms of r and t can be obtained by comparing the coefficient of 
(d9 2 + sin 2 6 dtp) in the two coordinate systems, leading to 

R(t,f)=r (27) 
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A second expression containing u, r and t can be found by making a transformation from LTB 
coordinates to find in null cone coordinates the term g 00 , which is equal to zero. A first order 
partial differential equation follows: 

Fin Fin 

= gf - ftf(t,r)]^ with «(t,0) = t (28) 

which, using the method of characteristics, can be written as: 

§ = -[i?,f(t,f)]. (29) 

Depending on the complexity of the R t ? term, either numerical or analytic methods can be used 
to solve this equation. 

Expressions containing partial derivatives for the other null cone terms can be obtained from 
covariant transformations. For h r and h u these are: 

Having values for h u and h r , the Bondi-Sachs coefficients can be obtained from: 

(3 = -\n\ h r \ and W — r(h u h r — 1) 

In order to write v\ in null coordinates, requires the transformation of the comoving velocity 
(v^ = (1,0,0,0) and v M = (—1,0,0,0)) into null coordinates. 

n du i dr dt dt 

«=^-,w=-^-,wo = —*- and «i = - — (32) 

at ot om ar 

The five point difference equations in |26| were used to determine the grid values of equations (|30|) . 
([5T|l and ([52^1 ; «i from (|3"2"]l together with p from ([231) are used as the input data on the initial null 
cone. 



C. Verification results 

A simple choice of bang function is implemented as a verification model: 

t B {r)=br, (33) 

with b being a constant. This simplifies equations (|23|) . (|24|) and (f25j) to: 

#(t,r) =r{t-br) 2/3 (34) 

l(-3i + 56r) 
^l*»H- 3 (< _ 6r) i/3 ^ 

^' r) = 2 ff (t-6r)(3t-56r))- ^ 

The effect of different values of b on the shape of the PNC is shown in figure @] where different 
models have been scaled and transposed relative to an observer located at the vertex of the PNC 
of a normalised EdS universe. The value b = is exactly the EdS model and b > shifts the 



11 

age of a universe to a younger age as r increases while b < provides the opposite effect where a 
universe is shifted to an older age as r increases. The latter case is particularly interesting since 
it provides a mechanism to mimic a cosmological constant [36| . In figure HI a reasonable match on 
the shape of the ACDM (0 A = 0.7) PNC is shown by a b = -0.5 bang time model before the PNC 
refocusses. This is however not an exact physical match but provides a useful verification model 
since the shape of the null cone is a critical aspect on the stability of a CIVP code. 




0.2 0.4 0.6 0.8 
Diameter distance (r) 



FIG. 4: The past null cone of different bang time models compared to the Einstein-de Sitter model and 
the ACDM model. 

The test cases presented in this section are the 6 = (EdS) case up to 0.82 u m ax and 0.82 r max 
and the b = —0.5 case up to 0.75 u max and 0.8 r max with u max and r max the time and distance 
of the AH. The EdS case provides a useful baseline since its solution is exactly known in null 
coordinates while the b = —0.5 case is representative of an LTB model with a ACDM-like PNC. 
It should be noted that the EdS model is not homeogeneous with respect to the null cone radial 
coordinate r and is therefore also representative of more general LTB models. The limits chosen are 
the maximum values where stable solutions were achieved. Extending these limits further causes 
the code breaks down rapidly, which is the expected behaviour close to the AH. 

The radial outer limits of the initial PNCs correspond to observations at z = 0.47 and z = 0.48 
in the EdS and LTB cases respectively. By reducing the extent of u limit, the radial extent 
can be increased and conversely reducing the radial limit allows the u limit to be increased, for 
instance values of 0.2 u max and 0.99 r max (z — 1.02) also produce stable results for the EdS 
model. The stability in terms of the radial and evolution limits is dependent on the input data and 
the simulations in section [V] were stable up to z « 1, which is a significant region for supernovae 
redshift-distance observations. In this section values were chosen that were useful for demonstrating 
accuracy and stability on normalised models. 



1. Physical variables 

Figure [5] shows the density (p) and covariant velocity (i>i) of the two models. The plots contain 
lines, which represent the exact values, and points representing the code values. The exact and 
code values are close enough not to show visible deviations. The density distributions of the two 
models on the oldest PNC are interesting to compare where the b = —0.5 model clearly displays 
the effect of the bang time surface. 
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FIG. 5: p and Vi against r on null cones in the past for the Einstein-de Sitter model (A) and (B) and the 
b — —0.5 bang time model (C) and (D). 



2. Local error propagation and convergence 



Figure [6] shows the radial distributions of the relative errors on the oldest null cones. As a 
measure of errors, the relative percentage of the difference between the exact and code of the 
density calculation has been chosen as a representative quantity. This choice has been based on 
different test cases where small errors in the other variables, v\, W and /3 have become pronounced 
in the density calculation. The EdS model having more accurate input values has a better local 
accuracy but the series-CIVP merger region is more visible in the r « region. 

In figure [3 errors are shown where grid resolution has been refined three times. It can be seen 
that errors in the region where the series solution is merged with the CIVP solution is significantly 
larger at low grid resolutions but converge rapidly to containable levels. The merging region be- 
tween the series and CIVP solution is particularly sensitive which suggests that further refinement 
of this region will be required especially if realistic observational data will be used as input in 
future research. 

In figure [8j the convergence against the radial grid size is shown. Here the slope of the line on 
a log scale is 1.9 and 1.8 for the EdS and b = —0.5 cases respectively, indicating that convergence 
approaches second order. The fact that the code breaks down after being close to second order 
accurate is as an accuracy consideration a useful feature since it doesn't produce false results; it 
either works accurately or provides no results. 
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FIG. 6: Spatial error vs. r for 200 r x 1000 u grid points for the Einstein-de Sitter model (A) and the 
b — —0.5 bang time model (B). 
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FIG. 7: Grid convergence: The error propagation for u x r grid resolutions of 1000 x 200, 500 x 100 and 
250 x 50 for the Einstein-de Sitter model (A) and the b — —0.5 bang time model (B). 
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FIG. 8: Grid convergence: The slope of the log(error max ) vs. log(Ar) line is 1.9 for a Ar/Au ss 3.1 grid 
for the Einstein-de Sitter model and 1.8 for a Ar/Au ss 3.2 grid for the b — —0.5 bang time model. 
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LTB VS. ACDM 



Overview 



Following from the Isotropic Observations Theorems in [9| (also see [30]), any reasonable set of 
spherical symmetrical observations can be accommodated by some distribution of inhomogeneities 
in an LTB model. Soon after the discovery of dark energy from Supernovae la observations in 
1998 [37, 38], using this principle, it was demonstrated by Celerier J36j], with others to follow 
(e.g. 39]), that the same observational effects can be caused by radial inhomogeneities without 
any need of a cosmological constant. The concept of inhomogeneities mimicking dark energy 
has been investigated many times since then, mostly as toy models demonstrating the ambiguity 
of observations commonly accepted in conventional cosmology. The main argument against LTB 
universes is that they place the observer in the centre of the Universe which, although not impossible 
is not a philosophically appealing idea. Even so, with a CIVP code at hand, an interesting numerical 
experiment is to test the historical evolution of an LTB model, with a zero cosmological constant, 
when observational data representing that of the ACDM model is interpreted as being caused by 
inhomogeneities . 

A simple simulation set up for this experiment is done by calculating the density p and covariant 
velocity v\ on the past null cone by transforming exact solutions for a flat ACDM model onto 
the past null cone. These are then used as input to the code and compared to the transformed 
model. In terms of the LTB metric (|21l) . using cosmological properties at the current epoch (to), 
the solution of the flat ACDM model is (see |40|): 



R(t,f) = S(t)f 






1/3 



sinh 



2/3 



(37) 



with the age of the Universe: 



h = 



sinh 



/ "A0 



in() 



1/2' 



(38) 



and the density distribution 



8^G S 3 ' 



(39) 



Here, fl m o and fl\o are the current density parameters for Baryonic matter and the cosmological 
constant respectively and H is the current Hubble constant. The values used to represent the 
actual Universe were: J7 m o = 0.3, Qao = 0.7 and Ho = 72 Mpc s _1 km -1 . An important issue to 
notice here is that p m in the ACDM model is determined by parameters related to the expansion 
and the density content and not by an independent measure of matter distribution such as number 
counts. The additional degree of freedom introduced in the LTB model is therefore not purely 
satisfied by an additional boundary condition which is a limitation that should be borne in mind 
when interpreting redshift dimming as an LTB model. 

Figure [S] shows the resulting LTB vs. ACDM evolutions back in time. While it might currently 
not be possible to distinguish these models from one another based on observations on the past null 
cone, in the past, these universes are distinctly different. In particular, the LTB universe seems 
to be heading towards a singularity much faster than the ACDM universe and would therefore be 
significantly younger, if the trend is to continue beyond the calculations. 
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FIG. 9: Density distribution (A) and covariant velocity (B) on past null cones at different proper times 
(u) 



B. LTB vs. ACDM matching 



A question arising from the results presented in figure [5] is: if the LTB null cones in the past 
would also represent ACDM flat space null cones? In this section an attempt is made to find a 
matching ACDM model for the to — 8Gyrs LTB null cone. The procedure followed here is to firstly 
find a null cone in the past of a selected ACDM model for which the density corresponds to the 
LTB density at r — 0. This gives a density curve which has the same starting point but the slope of 
the curve differs to that of the LTB model. The slope of the curve is then matched by adjusting the 
Hubble rate (-Ho) until it approximately corresponds to the LTB curve. The covariant velocities 
(vi as a function of r) are then compared to see, qualitatively, if the models can be regarded as 
the same. 

The models investigated are summarised in table U and detailed results are shown in the figures 
[TU]to[T3J From table U it can be seen that the matching instances are from universes with different 
ages and matching takes place at different times in their evolutions. From the detailed results, it 
becomes apparent that matching both the radial matter distribution and expansion through p and 
v\ is an unlikely proposition. Applying the matching procedure to p causes v\ to move away from 
the ACDM data while p moves away when v\ is matched. This then suggests that the past null 
cones do not represent a ACDM model. This is an illustration of similar conclusions found by Yoo, 
Kai & Nakao 0. 

An interesting correlation is that the values of Hq and t ma tch are very close for all the models 
and it might be possible that they are in fact the same. A formal correlation was however not done 
since the qualitative behaviour of the covariant velocities rules out a reasonable match in any of 
the models. The plots are presented in figures (J10[) to (| 13[) which reveal a surprising similarity in 
the vi plots. 



VI. CONCLUSION 



In this paper we demonstrated how to use the characteristic formalism of numerical relativity 
to investigate problems in observational cosmology. For this purpose, a CIVP code was developed 
and it was shown that the code is second order accurate and stable for selected LTB models in the 
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region before the PNC starts to refocus. By doing a numerical experiment with LTB vs. ACDM 
data, it was demonstrated that although the initial values of the two models can correspond on 
the current past null cone, the histories of the two models are distinctly different. The density of 
the LTB model rises significantly more quickly indicating a much younger universe, possibly too 
young. The result that past null cones evolved from an initial LTB null cone cannot be matched 
with a flat ACDM model has an important implication: While in our current epoch the LTB vs. 
ACDM ambiguity is difficult to disentangle, this is a feature of the Universe's current state and 
not its past. In other words, if the Universe is inhomogeneous without a cosmological constant, not 
only is the observer in a privileged position (near a central point), he also lives in a specific time 
where the Universe can appear to be either LTB or ACDM. It would be interesting to investigate 
extending this result to non-flat models. 
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